% calculate change in the weights given spike time lists which may 
% include negative spike times.  Assumes the lists are sorted.
%
% stdp_curve is the typical stdp curve.  It assumes that the middle is the
% time of the post-synaptic spike, and the left of middle is the LTP level
% for presynaptic spikes arriving before the post spike, and the right of
% middle is the LTD level for late arrivals.  Each entry is 1 ms delay.
%
% This code considers all spike interactions within the window, not just
% the nearest spikes.
% assumes w is normalized between 0 and 1.
function sd = stdp( prespiketimes_ms, postspiketimes_ms, delay_ms, stdp_curve, lrate, w)

% fix prespiketimes to deal with delay so I don't have to think about it
% again
prespiketimes_ms = prespiketimes_ms + delay_ms;


mid = floor(length(stdp_curve)/2)+1;


sd = 0;

% LTP
pre = 1;
for posttime = postspiketimes_ms(postspiketimes_ms>0)

    while pre <= length(prespiketimes_ms) & prespiketimes_ms(pre) < posttime
        pre = pre + 1;
    end
    % last pre spike is the closest one before the post spike.
    if( pre-1 > 0 )
        time_diff_index = mid + prespiketimes_ms(pre-1) -posttime ;
        if time_diff_index > 0 & time_diff_index <= length(stdp_curve)
            sd = sd + (1-w)*lrate* stdp_curve( time_diff_index );
%            disp([ stdp_curve(time_diff_index) time_diff_index]);
        end
    end
        
end

% LTD
post = 1;
for pretime = prespiketimes_ms(prespiketimes_ms>0)

    while post <= length(postspiketimes_ms) & postspiketimes_ms(post) <= pretime
        post=post+1;
    end
    % last pre spike is the closest one before the post spike.
    if( post-1 > 0 )
        time_diff_index = mid + pretime - postspiketimes_ms(post-1) ;
        if time_diff_index > 0 & time_diff_index <= length(stdp_curve)
            sd = sd + w* lrate* stdp_curve( time_diff_index );
%            disp([ stdp_curve(time_diff_index) time_diff_index]);
        end
    end
        
end




